MODULE pfmm_qlm_builder

   USE fmm_global_consts
   USE fmm_global_types
   USE pfmm_global_types
   USE fmm_stats
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: pfmm_get_raw_qlm,             &
!FIXME!!!!!!!!
             pfmm_deallocate_qlm

!FIXME!!!!!!
   ! Public variables
!   PUBLIC :: pfmm_system_size, pfmm_coord_shift

!   ! Coordinate translation required to ensure all primitive (x,y,z) positive
!   REAL(REALK), SAVE :: pfmm_coord_shift(3)
!   ! (minimum) dimension of a cube which encloses all moments
!   REAL(REALK), SAVE :: pfmm_system_size

   ! Number of electronic moments (LHS, RHS) and nuclear charges
   INTEGER(INTK), SAVE :: nTmoms, nSQmoms, nnuc
   ! Number of AO (contracted) basis functions
   INTEGER(INTK), SAVE :: nbas

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE pfmm_get_raw_qlm(scheme,dens,LHS,RHS,Tvec_map)

      IMPLICIT NONE
      TYPE(fmm_scheme),   INTENT(IN)  :: scheme
      REAL(REALK),        INTENT(IN)  :: dens(:,:,:)
      TYPE(raw_mm_data),  INTENT(OUT) :: LHS, RHS
      TYPE(T_vector_map), POINTER     :: Tvec_map(:)
      INTEGER(INTK) :: LMAX

      LMAX = scheme%raw_LMAX
      CALL pfmm_get_dims_from_file(scheme,LMAX)

      ! Allocate mulitpole moments arrays
      CALL pfmm_allocate_mm_arrays(LMAX,nTmoms,LHS)
      CALL pfmm_allocate_mm_arrays(LMAX,nSQmoms+nnuc,RHS)
      IF (ASSOCIATED(Tvec_map)) CALL fmm_quit('Tvec_map')
      ALLOCATE(Tvec_map(nTmoms))

      ! Read in moments from PMD1 code and store them in LHS and RHS
      CALL pfmm_read_nuclear_moms(RHS)
      CALL pfmm_read_electronic_moms(scheme,dens,LHS,RHS,Tvec_map)

      ! Stats and clean-up
      stat_raw_moms_LHS = SIZE(LHS%paras)
      stat_raw_moms_RHS = SIZE(RHS%paras)
      print *, "raw_moms_LHS =", SIZE(LHS%paras)
      print *, "raw_moms_RHS =", SIZE(RHS%paras)

   END SUBROUTINE pfmm_get_raw_qlm

!-------------------------------------------------------------------------------

   SUBROUTINE pfmm_get_dims_from_file(scheme,LMAX_in)

      IMPLICIT NONE
      TYPE(fmm_scheme), INTENT(IN) :: scheme
      INTEGER(INTK),    INTENT(IN) :: LMAX_in

      INTEGER(INTK) :: LMAX, ntot

      ! Read number of electronic moments
      nTmoms = 0
      nSQmoms = 0
      OPEN (UNIT=IO2,FILE='multipoles_header.md1',STATUS='OLD',   &
            ACTION='READ',FORM='UNFORMATTED')
      REWIND (IO2)
      READ (IO2) LMAX, nbas, nTmoms, nSQmoms
      CLOSE(UNIT=IO2,STATUS='KEEP')

      ! Read number of nuclear charges
      nnuc = 0
      OPEN (UNIT=IO2,FILE='nuclei.fmm',STATUS='OLD',   &
            ACTION='READ',FORM='FORMATTED')
      REWIND (IO2)
      READ (IO2,'(I8)') nnuc
      CLOSE(UNIT=IO2,STATUS='KEEP')

      ! Verifications
      WRITE(IO1,'(A,I8)') " Number of LHS electronic moments read =", nTmoms
      WRITE(IO1,'(A,I8)') " Number of RHS electronic moments read =", nSQmoms
      WRITE(IO1,'(A,I8)') " Number of nuclear charges read        =", nnuc
      ntot = nTmoms + nSQmoms + nnuc
      IF (ntot < 1) CALL fmm_quit('No moments generated!')
      IF (LMAX /= LMAX_in) THEN
         WRITE(IO1,*) LMAX, LMAX_in
         CALL fmm_quit('LMAX inconsistency in MM interface!')
      END IF

   END SUBROUTINE pfmm_get_dims_from_file

!-------------------------------------------------------------------------------

   SUBROUTINE pfmm_allocate_mm_arrays(LMAX,ndim,mm_arr)

      IMPLICIT NONE
      INTEGER(INTK),     INTENT(IN)  :: LMAX, ndim
      TYPE(raw_mm_data), INTENT(OUT) :: mm_arr
      INTEGER(INTK) :: i, foo, alloc_error

      NULLIFY (mm_arr%paras, mm_arr%dens, mm_arr%J_indices, mm_arr%batch_map)
      NULLIFY (mm_arr%qlm, mm_arr%qlm_W, mm_arr%qlm_T)

      ALLOCATE (mm_arr%paras(ndim))
      ALLOCATE (mm_arr%J_indices(ndim))
      ALLOCATE (mm_arr%dens(ndim))

      ! Initialize
      DO i = 1, ndim 
         mm_arr%paras(i)%cntr = zero
         mm_arr%paras(i)%ext = zero
         mm_arr%paras(i)%id = i
         mm_arr%paras(i)%batch = i
         mm_arr%paras(i)%map_up = 0
         mm_arr%paras(i)%box = 0 
         mm_arr%paras(i)%bra = 0
         mm_arr%paras(i)%box_cntr = zero
         mm_arr%J_indices(i)%i_indx = 0
         mm_arr%J_indices(i)%j_indx = 0
         mm_arr%dens(i) = zero
      END DO

      ! Full set of Qlm moments only relevant for electronic part
      foo = (LMAX+1)**2
      WRITE(IO1,*) 'mm_arrays: Attempting to allocate',  &
                     MAX(1,foo*ndim*8/1000000), 'MB of memory...'
      ALLOCATE (mm_arr%qlm(foo,ndim), STAT=alloc_error)
      IF (alloc_error /= 0) WRITE(IO1,*) '... Failed!'

      mm_arr%qlm(:,:) = zero  ! only non-zero written explicitly

   END SUBROUTINE pfmm_allocate_mm_arrays

!------------------------------------------------------------------------------

   SUBROUTINE pfmm_deallocate_mm_arrays(mm_data)

      IMPLICIT NONE
      TYPE(raw_mm_data), INTENT(OUT) :: mm_data
      INTEGER(INTK) :: i

      IF (ASSOCIATED(mm_data%paras)) DEALLOCATE (mm_data%paras)
      IF (ASSOCIATED(mm_data%J_indices)) DEALLOCATE (mm_data%J_indices)
      IF (ASSOCIATED(mm_data%dens)) DEALLOCATE (mm_data%dens)
      IF (ASSOCIATED(mm_data%qlm)) DEALLOCATE (mm_data%qlm)
      IF (ASSOCIATED(mm_data%qlm_T)) DEALLOCATE (mm_data%qlm_T)
      IF (ASSOCIATED(mm_data%qlm_W)) DEALLOCATE (mm_data%qlm_W)

      IF (ASSOCIATED(mm_data%batch_map)) THEN
         DO i = 1, SIZE(mm_data%batch_map)
            CALL free_batch_map(mm_data%batch_map(i)%head)
         END DO
      END IF
      IF (ASSOCIATED(mm_data%batch_map)) DEALLOCATE(mm_data%batch_map)

      NULLIFY (mm_data%paras, mm_data%dens, mm_data%batch_map)
      NULLIFY (mm_data%qlm, mm_data%qlm_W, mm_data%qlm_T)
      NULLIFY (mm_data%J_indices)

   CONTAINS 

      RECURSIVE SUBROUTINE free_batch_map(node)

         IMPLICIT NONE
         TYPE(id_node), POINTER :: node
         IF (ASSOCIATED(node%next)) THEN
            CALL free_batch_map(node%next)
         END IF
         DEALLOCATE(node)
         NULLIFY(node)

      END SUBROUTINE free_batch_map

   END SUBROUTINE pfmm_deallocate_mm_arrays

!-------------------------------------------------------------------------------

   SUBROUTINE pfmm_deallocate_qlm(LHS,RHS)

      IMPLICIT NONE
      TYPE(raw_mm_data), INTENT(OUT) :: LHS, RHS

      CALL pfmm_deallocate_mm_arrays(LHS)
      CALL pfmm_deallocate_mm_arrays(RHS)

   END SUBROUTINE pfmm_deallocate_qlm

!-------------------------------------------------------------------------------

   SUBROUTINE pfmm_read_nuclear_moms(mm_arr)

      IMPLICIT NONE
      TYPE(raw_mm_data), INTENT(OUT) :: mm_arr
      REAL(REALK)   :: PX,PY,PZ
      INTEGER(INTK) :: Q, i

      ! Next read nuclei data: charge and location
      ! Positive charges give a positive potential contribution
      IF (nnuc == 0) RETURN
      OPEN (UNIT=IO2,FILE='nuclei.fmm',STATUS='OLD',  &
            ACTION='READ',FORM='FORMATTED')
      REWIND (IO2)
      READ (IO2,'(I8)') i

      WRITE (IO1,*) "Reading in point charges:"
      DO i = 1, nnuc 
         READ (IO2,'(F25.15)') PX
         READ (IO2,'(F25.15)') PY
         READ (IO2,'(F25.15)') PZ
         READ (IO2,'(I5)') Q
         WRITE (IO1,'(I4,X,3F8.4,I4)') I, PX,PY,PZ, Q
         mm_arr%qlm(1,i) = DBLE(Q)
         mm_arr%paras(i)%cntr = (/ PX, PY, PZ /)
      END DO

      mm_arr%dens(:) = one

      CLOSE(UNIT=IO2,STATUS='KEEP')
!FIXME
!      CALL pfmm_get_system_size_and_shift(mm_nuc%paras)

   END SUBROUTINE pfmm_read_nuclear_moms

!-------------------------------------------------------------------------------
! Read in multipole moment data from interface file.
! In all this MM code we assume the order of moments is:
!   (0),(-1,0,1),(-2,-1,0,1,2)... wrt (L,M)
!
! Also we assume in this routine that the qlm are written in batches
! The first batch = (mu^0,nu^T).
! The second      = (lambda^S,sigma^S+Q1) for all S.
! The third       = (lambda^S,sigma^S+Q2) for all S.
! etc
!
! End-of-batch points are denoted by special L-values as commented below

   SUBROUTINE pfmm_read_electronic_moms(scheme,dens,LHS,RHS,Tvec_map)

      IMPLICIT NONE
      TYPE(fmm_scheme),  INTENT(IN)  :: scheme
      REAL(REALK),       INTENT(IN)  :: dens(:,:,:)
      TYPE(raw_mm_data), INTENT(OUT) :: LHS, RHS
      TYPE(T_vector_map), POINTER    :: Tvec_map(:)
      REAL(REALK)   :: PX,PY,PZ, SPH
      INTEGER(INTK) :: I,L,M, A,B, LM, nmom, Tindex, Qindex
      INTEGER(INTK) :: ivec(3)

      ! Read electronic multipole moments into core
      IF (nTmoms + nSQmoms == 0) RETURN
      OPEN (UNIT=IO2,FILE='multipoles.md1',STATUS='OLD',  &
            ACTION='READ',FORM='UNFORMATTED')
      REWIND (IO2)


      ! Dummy first line for T-batch
      WRITE (IO1,*) "Reading in (mu,nu) moms:"
      READ (IO2) I,L,M, A,B, PX,PY,PZ, SPH
      IF (L /= -2) CALL fmm_quit('no PFMM T-batch!')
      nmom = 0
      readloop1: DO
         ! Read T-vector
         READ (IO2) Tindex, ivec(1), ivec(2), ivec(3)
print *, "T-vec(:) =", Tindex, ivec(:)
         readloop2: DO
            READ (IO2) I,L,M, A,B, PX,PY,PZ, SPH
            ! End of all T-batches marked by record with L=-22
            IF (L == -22) EXIT readloop1
            ! End of this T-batch marked by record with L=-2
            IF (L == -2) EXIT readloop2
            nmom = nmom+1
            CALL pfmm_store_mm(Tindex,nmom,L,M,A,B,PX,PY,PZ,SPH,dens,LHS)
            Tvec_map(nmom)%Tvec(:) = ivec(:)
            Tvec_map(nmom)%Tindex = Tindex
         END DO readloop2
      END DO readloop1


      ! Dummy first line for SQ-batches
      WRITE (IO1,*) "Reading in (lambda,sigma) moms:"
      READ (IO2) I,L,M, A,B, PX,PY,PZ, SPH
      IF (L /= -3) CALL fmm_quit('no PFMM SQ-batches!')
      nmom = 0
      readloop3: DO
         ! Read Q-vector
         READ (IO2) Qindex, ivec(1), ivec(2), ivec(3)
print *, "Q-vec(:) =", Qindex, ivec(:)
         readloop4: DO
            READ (IO2) I,L,M, A,B, PX,PY,PZ, SPH
            ! End of one SQ-batch marked by record with L=-3
            IF (L == -3) EXIT readloop4
            ! End of all SQ-batches marked by record with L=-33
            IF (L == -33) EXIT readloop3
            nmom = nmom+1
            CALL pfmm_store_mm(Qindex,nmom+nnuc,L,M,A,B,PX,PY,PZ,SPH,dens,RHS)
         END DO readloop4
      END DO readloop3


      CLOSE(UNIT=IO2,STATUS='KEEP')
!      CALL pfmm_get_system_size_and_shift(mm_elec%paras)

   END SUBROUTINE pfmm_read_electronic_moms

!-------------------------------------------------------------------------------

   SUBROUTINE pfmm_store_mm(Dx,I,L,M,A,B,PX,PY,PZ,SPH,dens,mm_arr)

      IMPLICIT NONE
      INTEGER(INTK),     INTENT(IN)  :: Dx
      INTEGER(INTK),     INTENT(IN)  :: I,L,M, A,B
      REAL(REALK),       INTENT(IN)  :: PX,PY,PZ, SPH
      REAL(REALK),       INTENT(IN)  :: dens(:,:,:)
      TYPE(raw_mm_data), INTENT(OUT) :: mm_arr

      INTEGER(INTK) :: LM

      IF (ABS(SPH) > 1.0d-12)  &
        WRITE (IO1,'(5I4,X,3F8.4,2E13.4)') I,L,M,A,B,PX,PY,PZ,SPH,dens(A,B,Dx)

      IF (A > nbas) CALL fmm_quit('pfmm_store_mm error 0')
      IF (B > nbas) CALL fmm_quit('pfmm_store_mm error 00')
      IF (I < LBOUND(mm_arr%qlm,2)) CALL fmm_quit('pfmm_store_mm error 11')
print *, I, UBOUND(mm_arr%qlm,2)
      IF (I > UBOUND(mm_arr%qlm,2)) CALL fmm_quit('pfmm_store_mm error 111')
      IF (((L+1)**2) > SIZE(mm_arr%qlm,1)) CALL fmm_quit('interface 156')

      ! Indices to map moments to orbitals in J-matrix
      mm_arr%J_indices(I)%i_indx = A
      mm_arr%J_indices(I)%j_indx = B

      ! Density matrix factor
      IF (Dx > UBOUND(dens,3)) CALL fmm_quit('interface error Dens')
      mm_arr%dens(I) = -dens(A,B,Dx)

      ! Multipole expansion centre
      mm_arr%paras(I)%cntr = (/ PX, PY, PZ /)

      ! See defn of 'extent' in p.424 MEST Helgaker et al
      mm_arr%paras(I)%ext = 0

      ! Components (l,m) of multipole without density factorized in 
      LM = L*(L+1) +M +1 
      mm_arr%qlm(LM,I) = SPH

   END SUBROUTINE pfmm_store_mm

!-------------------------------------------------------------------------------
!
!   SUBROUTINE pfmm_get_system_size_and_shift(paras)
!
!      IMPLICIT NONE
!      TYPE(raw_mm_paras), INTENT(INOUT) :: paras(:)
!      REAL(REALK)   :: sys_min(3), sys_max(3)
!      INTEGER(INTK) :: i
!
!      sys_min = paras(1)%cntr
!      sys_max = paras(1)%cntr
!      DO i = 1, SIZE(paras)
!         sys_min(:) = MIN(sys_min(:),paras(i)%cntr(:))
!         sys_max(:) = MAX(sys_max(:),paras(i)%cntr(:))
!      END DO
!      pfmm_system_size = MAXVAL(sys_max - sys_min)
!!      IF (pfmm_system_size == zero) CALL fmm_quit('zero system size!')
!
!      pfmm_coord_shift = sys_min
!!      DO i = 1, SIZE(paras)
!!         paras(i)%cntr(:) = paras(i)%cntr(:) - sys_min(:)
!!      END DO
!
!   END SUBROUTINE pfmm_get_system_size_and_shift
!
!-------------------------------------------------------------------------------
!
!   SUBROUTINE pfmm_assign_LHS_RHS_data(LMAX,mm_elec,mm_nuc,LHS,RHS)
!
!      IMPLICIT NONE
!      INTEGER(INTK),     INTENT(IN)  :: LMAX
!      TYPE(raw_mm_data), INTENT(IN)  :: mm_elec
!      TYPE(raw_mm_data), INTENT(IN)  :: mm_nuc
!      TYPE(raw_mm_data), INTENT(OUT) :: LHS, RHS
!
!      INTEGER(INTK) :: i
!
!      ! Assign LHS data (mu,nu^T) AO pairs
!      CALL pfmm_allocate_mm_arrays(LMAX,nTmoms,LHS)
!      ! Only electronic
!      LHS%qlm(:,:)     = mm_elec%qlm(:,1:nTmoms)
!      LHS%dens(:)      = zero
!      LHS%paras(:)     = mm_elec%paras(1:nTmoms)
!      LHS%J_indices(:) = mm_elec%J_indices(1:nTmoms)
!      ! Re-scaled moments
!      NULLIFY(LHS%qlm_T)
!      NULLIFY(LHS%qlm_W)
!
!      ! Assign RHS data (lambda^S,sigma^S+Q) AO pairs + nuclei
!      CALL pfmm_allocate_mm_arrays(LMAX,nSQmoms+nnuc,RHS)
!      ! Electronic first
!      LHS%qlm(:,1:nSQmoms)     = mm_elec%qlm(:,1:nSQmoms)
!      LHS%dens(1:nSQmoms)      = mm_elec%dens(1:nSQmoms)
!      LHS%paras(1:nSQmoms)     = mm_elec%paras(1:nSQmoms)
!      LHS%J_indices(1:nSQmoms) = mm_elec%J_indices(1:nSQmoms)
!      ! Nuclear next
!      LHS%qlm(:,nSQmoms+1:)     = zero
!      LHS%qlm(1,nSQmoms+1:)     = mm_nuc%qlm(1,:)
!      LHS%dens(nSQmoms+1:)      = mm_nuc%dens(:)
!      LHS%paras(nSQmoms+1:)     = mm_nuc%paras(:)
!      LHS%J_indices(nSQmoms+1:) = mm_nuc%J_indices(:)
!      ! Re-scaled moments
!      NULLIFY(RHS%qlm_T)
!      NULLIFY(RHS%qlm_W)
!
!      ! Initialise parameter:moments mapping, 
!      ! and batch numbers for unique centres
!      DO i = 1, SIZE(LHS%paras)
!         LHS%paras(i)%id = i
!         LHS%paras(i)%batch = i
!      END DO
!      DO i = 1, SIZE(RHS%paras)
!         RHS%paras(i)%id = i
!         RHS%paras(i)%batch = i
!      END DO
!
!   END SUBROUTINE pfmm_assign_LHS_RHS_data
!
!-------------------------------------------------------------------------------

END MODULE pfmm_qlm_builder

